//VerificationSim.h
#include <vector>
#include <iostream>
#include <math.h>
#include <ctime>
#include <string>
#include <cstdlib>
using namespace std;
//std::vector< std::vector<double> > A;

extern const double inverse_sqrt2 = 0.707106781;
extern bool trace;


class Matrix{
	
		vector< vector<double> > contents; 
	
	public:
		int size; //size of one side of the matrix. Number of elements is n^2
	
	//CONSTRUCTORS
		//default constructor
		Matrix(){
		};
		
		//copy constructor
		Matrix(const Matrix& a){
			size = a.size;
		};
		
		//default constructor for size n
		Matrix(int n){
			size = n;
			contents = vector< vector<double> > (n, vector<double>(n));
		};

		Matrix(int n, bool v){
			double input;
			size = n;
			contents = vector< vector<double> > (n, vector<double>(n));
			for(int a = 0; a < n; a++){
				for(int b = 0; b < n; b++){
					cin >> input;
					contents[a][b] = input;
				}
			}
		}

		
	//FUNCTIONS
		
		//Returns the value at the matrix location
		double Get(int row, int col){
			return contents[row][col];
		};
		
		//Sets a specific value in the matrix
		void Set(int row, int col, double num){
			contents[row][col] = num;
		};
		
		//Returns the contents of the matrix, without the matrix
		vector< vector<double> > getContents(){
			return contents;
		};
		
		//Print function, to output matrices. 
			//will not necessarily format well if numbers with varying number of digits are used.
		void Print(){
			cout << endl << "[";
				for(int a = 0; a < size;){
					for(int b = 0; b < size;){
						cout << Get(a,b);
						b++;
						if(b < size)
						cout << ", ";
					}
					a++;
					if(a < size)
					cout << endl << " ";
				}
			cout << "]" << endl;
		};
		
		//input function for construction, used with copy constructor
		static Matrix InputMatrix(int n){
			Matrix A(n);
			double input;
			for(int a = 0; a < n; a++){
				for(int b = 0; b < n; b++){
					cin >> input;
					A.Set(a,b,input);
				}
			}
			return A;
		};

		static Matrix SetPauliZ(){
			Matrix temp = Matrix(2);
			temp.Set(0,0,1);
			temp.Set(0,1,0);
			temp.Set(1,0,0);
			temp.Set(1,1,-1);
			return temp;
		};
};


extern Matrix* A;
extern Matrix* B;
extern Matrix* C;

class Qbit{
public:
	double re0;
	double im0;
	
	double re1;
	double im1;
	
//CONSTRUCTORS
	//default constructor
	Qbit(){
		re0 = 1;
		im0 = 0;
		re1 = 0;
		im1 = 0;
	};
	
	//copy constructor
	Qbit(const Qbit& a){
		re0 = a.re0;
		im0 = a.im0;
		re1 = a.re1;
		im1 = a.im1;
	};
	
	//Normal constructor
	Qbit(double _re0, double _im0, double _re1, double _im1){
		re0 = _re0;
		im0 = _im0;
		re1 = _re1;
		im1 = _im1;
	};
	
//FUNCTIONS
	//Normalizes the state, setting it to the normalized version
	void Normalize(){
		double total = pow(re0,2) - pow(im0,2) + pow(re1,2) - pow(im1,2);
		re0 = re0/total;
		im0 = im0/total;
		re1 = re1/total;
		im1 = im1/total;
		//return Qbit(re0,im0,re1,im1);
		
	}; //NOTE: I can change this to, instead, return the normalized version. If the original state is still needed.
	
	//Measures the qubit, and sets it to strictly |0> or |1>
	void Measure(){
		Normalize();
		double prob0 = pow(re0,2) - pow(im0,2);
		double prob1 = pow(re1,2) - pow(im1,2);
		if(prob0 > prob1){
			re0 = 1;
			im0 = 0;
			re1 = 0;
			im1 = 0;
		}
		else{
			re0 = 0;
			im0 = 0;
			re1 = 1;
			im1 = 0;
		}
	};
	
	void ApplyHgate(){
		double _re0 = inverse_sqrt2 * re0 + inverse_sqrt2 * re1;
		double _im0 = inverse_sqrt2 * im0 + inverse_sqrt2 * im1;
		double _re1 = inverse_sqrt2 * re0 - inverse_sqrt2 * re1;
		double _im1 = inverse_sqrt2 * im0 - inverse_sqrt2 * im1;
		re0 = _re0;
		im0 = _im0;
		re1 = _re1;
		im1 = _im1;
	};
	
	void ApplyGate(Matrix * gate){
		double _re0 = gate->Get(0,0) * re0 + gate->Get(0,1) * re1;
		double _im0 = gate->Get(0,0) * im0 + gate->Get(0,1) * im1;
		double _re1 = gate->Get(1,0) * re0 + gate->Get(1,1) * re1;
		double _im1 = gate->Get(1,0) * im0 + gate->Get(1,1) * im1;
		re0 = _re0;
		im0 = _im0;
		re1 = _re1;
		im1 = _im1;
	};
	
	void Print(){
		if(re0 > 0 && im0 > 0){
			if(re0 > 0){
				cout << re0;
				if(im0 > 0){
					cout << "+";
				}
			}
			if(im0 > 0){
				cout << im0 << "i";
			}
			cout << "|0>";
		}
		if(re1 > 0 && im1 > 0){
			if(re1 > 0){
				cout << re1;
				if(im1 > 0){
					cout << "+";
				}
			}
			if(im1 > 0){
				cout << im1 << "i";
			}
			cout << "|1>";
		}
		
	};
};

class Qstate{
	vector<Qbit> bits;
	
	Qstate(){
	};
	Qstate(Qbit temp){
		bits.insert(bits.end(), temp);
	};
	Qstate(vector<Qbit> a){
		bits = a;
	};
	
	Qbit Access(int location){
		return bits[location];
	};
};

Qbit initialization(Matrix* A, Matrix* B, Matrix* C, int n, int l, int k);


